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Abstract 

I report on a numerical program for the evolution of parton distributions. The program uses 
the Mellin-transform method with an optimized contour. Due to this optimized contour the 
program needs only a few evaluations of the integrand and is therefore extremely fast. In 
addition, the program can also be used to reproduce the results of the x-space method. 



PROGRAM SUMMARY 

Title of program: partonevolution 
Version: 1.0 
Catalogue number: 



Program obtained from: [http : / /www . f is .unipr . it/~stefanw/partonevolution 

E-mail: stefanw@fis.unipr.it 

License: GNU Public License 

Computers: all 

Operating system: all 

Program language: C++ 

Memory required to execute: negligible (1 MB) 
Other programs called: none 
External files needed: none 

Keywords: Parton distributions, next-to-leading order evolution. 

Nature of the physical problem: Parameterizations of parton distributions are usually given 
at a fixed input scale Qq. To be used in perturbative calculations, they have to be evolved 
numerically to the desired scale Q. 

Method of solution: Inverse Mellin transform method. 

Restrictions on complexity of the problem: None. 

Typical running time: lO^^s on a standard PC for an accuracy of 0.02%, see also sect. ^ 
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LONG WRITE-UP 



1 Introduction 

Experiments at hadron colliders require apart from the hard scattering amplitudes, which can 
be calculated in perturbation theory, also the knowledge of the parton distributions functions 
(pdf's). Although these parton distribution functions parameterize the non-perturbative infor- 
mation, they depend on a factorization scale and the dependence on that scale can be calculated 
within perturbation theory. Usually parameterizations of parton distributions functions are given 
at a low input scale [|T]] - ^ and then evolved to the desired scale. At present the evolution ker- 
nel is known completely to next-to-leading order (NLO) 0] - . The evolution has to be done 
numerically. There are two methods available: The first one consists in integrating numerically 
the evolution equations in x and in and is called the x-space method. The second one first 
performs a Mellin-transform on the variable x [[T^, [TT]] . The evolution equations then factorize 
and can be solved analytically. Using this Mellin-transform method (or N-space method), leaves 
us with the task to evaluate the inverse Mellin transform numerically. In this paper I describe 
a numerical program, which uses an optimized integration contour in the complex N-plane in 
order to do the inverse Mellin transform. This method has been proposed by Kosower [fT^]. 
Using this method gives a fast and accurate program for the evolution of parton distributions. 
Although a variety of numerical evolution programs already exist (mainly within the x-space 
method) [[Dp - p^], they usually cannot compete in computation speed. The required CPU time 
becomes an issue if one tries to extend the analysis of experiments at hadron colliders in two 
important directions: 

• Parton distributions with error-bars [ ^3] , ^ . The Tevatron data on the measurement of 



the one jet inclusive transverse energy has shown p^ ] that the method of obtaining par- 
ton distributions from "global fits" has reached its limits. Giele, Keller and Kosower 



p3| , p4| ] have proposed a method to include uncertainties into the parton distributions. 
This method involves a functional integration over a set of parton distributions. With 
a set of pdf's containing 1000 or more parton distributions, the traditional approach of 
a pre-calculated grid containing the evolved parton distributions becomes inappropriate. 
Instead, one would need a program which allows to calculate the evolution on the fly. The 
program presented here allows one to do just that. 

• Extension to NNLO. The calculation of the three-loop anomalous dimensions is currently 
under way p^ , [27] ]. It is likely that the result will be rather lengthy. The computational 
cost of evaluating the anomalous dimensions will therefore be rather high, and one would 
like to minimize the number of function calls to the anomalous dimensions. The opti- 
mized contour employed in this program minimizes this number. Once the three-loop 
anomalous dimensions are known, the program can be extended in a straightforward way 
to include NNLO corrections. 

The program is optimized for the N-space method, where one needs typically only four or five 
evaluations of the anomalous dimensions to obtain the evolution of a parton distribution with 
an accuracy of 2 ■ 10^. In addition, I also include an option to simulate the x-space evolution 
programs, x-space programs use a different truncation prescription and the results might differ 
numerically from the N-space approach. It should be noted, that both approaches are valid to 
next-to-leading order, the difference is formally of higher order. The x-space option is included 
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to allow an easy comparison with those programs, however it is not optimized. In technical 
terms, since no closed formula for the evolution operator in the singlet sector within the x-space 
approach is known, the evolution operator is approximated by a sequence of evolutions over a 
small interval. This slows down the program considerably. The N-space evolution method is 
therefore the recommended one, where the full performance of the program is obtained. 

This paper is organized as follows. In section |^ I recall the basic formulae for the evolution 
of parton densities and sketch how to choose an optimized contour using the Mellin transform 
method. Section ^ gives an overview of the design of the program, while section ^ is of a 
more practical nature and describes how to install and use the program. In section ^ I com- 
pare the program with the reference results of Bliimlein et al. and address issues like 
performance. A summary is provided in section ^ The appendix collects the formulae for the 
evolution operators. 



2 Theoretical Background 

The evolution equations for the parton distributions f{x, Q^) of the proton read 

2'^%^ = P{x,Q')®f{x,Q'), (1) 

where x stands for the nucleon's momentum fraction carried by the parton, P{x,Q^) is the 
Altarelli-Parisi evolution kernel, and ® denotes the convolution 

1 1 

A{x)®B{x) = j dy j dzb{x-yz)A{y)B{z). (2) 



Eq. (HI) is understood to represent the quark non- singlet evolution as well as the coupled quark 
singlet and gluon evolution. In the latter case / is a two-component vector and P{x, Q^) a 
two-by-two matrix. Eq. ([1]) can be factorized by taking Mellin moments 

1 

A~ = j dx?^-^A{x) (3) 



and one obtains 

2'^^^ = P'{Q^)-f{Q^). (4) 

Changing variables by using = a, (2^) /47l instead of In 2^ as evolution variable one obtains 



where 



or 
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|3(a.) = 2^3^ (6) 



is the beta function for the strong coupling. Up to now, the evolution kernel is known completely 
only to next-to-leading order: 



P^Q") = asPl + a]Pl + 0{a]). (7) 
To this order, the beta function is given by 

|3(a,) = -^Qa]-^ia] + 0{a1). (8) 
Inserting eq. and eq. (||) into eq. @ and expanding the r.h.s in Us one obtains 



das Po«.5 



f{Q^). (9) 



Solutions of eq. (|^) are referred to as "N-space" solutions, since this is the standard approach 
within the N-space method. One the other hand, one might defer from inverting the power series 
in the denominator on the r.h.s of eq. and to solve eq. (g) directly. To NLO, this amounts to 
solving 

M = -J^±fi^/=(e). (10) 

das as{^o + ^ias) 

Solutions of eq. ( [T0| ) are referred to as "x-space" solutions, again since this corresponds to the 
usual procedure within the x-space method. Eq. (g) and eq. (|TD|) differ by terms which are 
formally of higher order in as, and solutions to both equations are therefore considered accurate 
to next-to-leading order, although they might differ numerically. 

In addition there are different common practices on how to calculate at the scale Q^. At 
NLO as{Q^) can be found from the solution of the equation 

l-eil„(p, + 5£)_p„i = 0, (11) 

as po V ^.^y 



where L = In-^ — . This equation has to be solved numerically, and solutions of eq. ( |1 ID are 

^QCD 

referred to as "exact" NLO solutions. On the other hand it is common practice to view as as a 
power series in 1/L and to use 

1 / pi lnL\ 
= ^ (12) 



to calculate as at the scale Q^. Solutions of eq. ([T2|) are referred to as "truncated in 1 /L" solu- 
tions. Again, the numerical values might differ for solutions of eq. ([Tl|) and eq. ([T2|). It should 
be mentioned that the integration constant Aqcd appearing in eq. ([iT]) is not identical to the 
integration constant of the same name in eq. (jT^, e.g. fixing as at a scale like m| and solving 
eq. ( [TT| ) and eq. ( [T2| ) for Aqcd will yield two different numerical values. Each value has to be 
used in conjunction with the equation from which it was obtained. 

Eq. and eq. ( PUj ) are formally solved by the introduction of an evolution operator E^: 

f{Q^) = E'iasiQ^)MQl))ml) (13) 



5 



Furthermore, there is yet another truncation method used in the literature: Within this method 
one does not perform the change of the evolution variable from to a.,, but one substitutes the 



truncated solution eq. (12) for cis into eq. (H). This amounts to solving 



,23/=(e') / w, Pi inL\ ^ , 1 




I will refer to solutions of eq. (jl^) as "x-space and truncated in 1/L" solutions. Again, this 
equation is formally solved by an evolution operator, which now depends on and Qq instead 
ofusiQ^) and asiQl): 

m^) = E\Q\Qi)nQi) (15) 

Formulae for the evolution operators corresponding to the various truncation prescriptions are 
collected in the appendix. The evolved parton distributions in x-space are then obtained by an 
inverse Mellin transformation: 

c 

Here, the contour C runs to the right of all singularities of the integrand. 

The parton distributions are usually parametrized at the input scale in a form like 

xfMl) = £a,x«'(1-x)P' (17) 



with Mellin transform 



nal) = £A,5(2 + a,-l,l + |3,), (18) 



where B{x,y) is Euler's beta function. Dropping from now on the arguments and our task 
is to evaluate the integral 

/ = Re - fdzE^F{z), 
ni J 

Cs 

Fiz)=x~' £A,-5(z + a, - 1, 1 + |3,-), (19) 



where complex conjugation has been used to replace the contour C by Cs, starting at the real 
axis, right to the right-most pole and running upwards to infinity. In the singlet case F{z) and 
A, are vectors and is a matrix. The most efficient way to solve this problem is to choose a 
contour in such a way that the integrand can very well be approximated by some set of orthog- 
onal polynomials. The method is due to Kosower [[T^]. It uses a parabolic contour, completely 
determined by the parton distribution at the original scale Qq and evaluates the integral by a 
Gauss-Laguerre quadrature formula. With this method a few (e.g. four or five) evaluations of 
the integrand are sufficient to evaluate the integral to very high precision (e.g. to an accuracy of 
0.02 %). 

I shortly review this method for the non-singlet case. The contour of integration is given by 

z{u) ^ Z0 + iC2\/u + ^C2C3U, M = O...00, (20) 
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where zo is the minima of the function F{z) on the real axis right to the rightmost pole of F{z) ■ 
The parameters C2 and C3 are given by 



C2 



' 2F{zo) 
F"izo) ' 



C3 



F"'{zo) 
3F"{zoy 



(21) 



F'{z), F"{z) and F"\z) are the first three derivatives of F{z) with respect to z- For the parame- 
terization eq. (jl7|) they are given by 

F" = 1^{G': + G?)F,, 

i 

F'" = £(Gr + 3GfG;. + Gf)i^-, (22) 



where Fi=x ^AiB{z + a, - 1 , 1 + P/) and 



G, 
G'l 



ln(jc-^5(z + a,-l,l + (3,)), 
- Inx + + a, - 1 ) - + a/ + p/ 
V (z + a,- - 1 ) - V (z + a,- + p,-) , 
\|/'(z + a/- 1)- V(z + a/ + p/). 



(23) 



Here, \|/, \|/' and \[/" are polygamma functions. With this parameterization one obtains for the 
integral eq. ([TP|): 



C2 
271 



e "Re 



e"(l-/C2C3v^)£^(")F(z(i<)) 



(24) 



The function the weight function of the Laguerre polynomials ^^'^{x) and the 

integral can be approximated by a Gauss-Laguerre quadrature formula: 



^ £ Re (1 - ic2C,^) E'^"J^F {ziuj)) 



(25) 



— 1/2 

Here, are the zeros of (x) and the weights are given by 

r(n + i) 



(26) 



In the singlet case one has a two-component vector given by the quark singlet combination 
and the gluon distribution G^. There will be strong mixing already for moderate evolution 
in Q^, since and G^ are not eigenfunctions of the evolution kernel. One might be tempted 
to work in a basis of eigenfunctions of the leading-order evolution kernel and to extend the 
method discussed above to these eigenfunctions, eventually choosing two different contours for 
the two eigenfunctions. However this approach faces two problems: First of all, one of the two 
eigenfunctions usually does not have a minimum on the real axis on the right of the rightmost 



7 



pole. Secondly, by using different contours for the two components one is forced to keep track 
of branch cut crossings, a task one would like to avoid. It is therefore better to use a single 
contour for both components, although it might not be the optimal one for both components. 
In the singlet sector the contour is chosen as follows (this differs from the original proposition 
of Kosower): zo is chosen as the minimum of the sum of the E^- and the G^-component. I then 
diagonalize the leading-order evolution kernel at zo- 



{zo)-'k+{zo) -yO,gq{zo) 




(27) 



where the eigenvalues X± are given by 



(zo) + f^g' (zo) ± \l ( f^g' (zo) - fqcl (zo) ) + 4fqg' iZ0)f^q^ (zo) 



JO) I 



,(0) 



(0), 



(28) 



Note that any z-dependence in eq. ( p7| ) is only due to E' or G', zo appearing in the anomalous 



dimensions is a fixed parameter. Strictly speaking the eigenvectors of eq. (^7|) are only eigen- 
vectors at the point z = zo of the leading-order evolution kernel. I then take e\ to determine C2 
and C3 and integrate both components with this contour. To compensate for the fact that zo is 
not necessarily the minimum of e\, the formula eq. ( Plj ) is replaced by 



C2 



2F{zq)F"{zq) 



F"{zor-F'{zo)F"'izoy 



(29) 



which takes into account corrections due to the fact, that the contour does not start from a 
minimum of F{z) (e.g. F'{zo) ^ 0). The reasons why this works is that in the cases of interest e\ 
is roughly the sum of and G^ and so zo is approximately optimal for e\ . The other component 
^2 is roughly the difference and numerically therefore smaller. 



3 Design of the program 

The program is written in C-i-i-. All classes and functions are defined in the namespace pdf . 
The parameterization of a parton distribution at the initial scale Qq is represented by the class 
partondistribution. A partondistribution is constructed as follows: 

// u-valence CTEQ 4M 
int n = 2; 
double QO = 1 . 6; 

double A_u[2] = { 1.344, 1.344*6.402}; 

double alpha_u[2] = { 0.501, 0.501+0.873 } ; 
double beta_u[2] = { 3.689, 3.689 } ; 
int eta = -1; 

partondistribution uvalence = 

partondistribution (n, QO, A_u, alpha_u, beta_u, eta) ; 

This corresponds to the parameterization 

xu.ixMl) = 1.344;c0-50i(l-x)3-689(i+6.402/-«^3) ^^^^ 
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at the initial scale Qo= 1.6 GeV. Valence like distributions need to have T] = — 1 explicitly in the 
constructor. For distributions which correspond to T] = 1 (e.g. distributions in the singlet sector 
or non-singlet non-valence distributions like (m + m) — (d + d)), the rj-parameter can be dropped. 

This partondistribution can be evolved to a scale Qby calling the function evolve_nonsinglet: 

// evolution to 10 GeV for x = 0.01 
double Q = 10.0; 
double X = 0.01; 

double u_v = evolve_nonsinglet (uvalence, Q, x) ; 

The function evolve_nonsinglet is just an entry point and the job is passed to the function 
evolve_nonsinglet_5, which performs a Gauss-Laguerre quadrature with 5 points. In addi- 
tion there are functions evolve_nonsinglet_3, evolve_nonsinglet_10, 
evolve_nonsinglet_2 and evolve_nonsinglet_30, corresponding to quadrature formulae 
with 3, 10, 20 and 30 points. These functions can be selected by passing an additional argument 
to evolve_nonsinglet, e.g. 

// evolution using a 10-point quadrature formula 
double u_v = evolve_nonsinglet (uvalence, Q, x, 10) ; 

will call the function evolve_nonsinglet_10. 

In the singlet sector evolution is done with the help of the function evolve_singlet: 
// singlet evolution 

vector_d singlet = evolve_singlet ( Sigma, Gluon, Q, x) ; 

double s = singlet [0]; 
double g = singlet [1]; 

Sigma and Gluon are of type partondistribution and represent the parameterization of the 
quark singlet combination and the gluon distribution, respectively. The evolution function re- 
turns a two-component vector, with the evolved quark singlet distribution as first component 
and the evolved gluon distribution as second component. In C and C-i-i- subscription of arrays 
starts at zero. As in the non-singlet case the function evolve_singlet is just an entry point, 
calling evolve_singlet_5. Again, there are functions using quadrature formulae with 3, 5, 
10, 20 and 30 points. 

In addition there are the classes evolutionkernel and alpha_strong. The class 
evolutionkernel has a static data member method, which allows one to select the truncation 
method, e.g. 

// select truncation method for evolution kernel 
evolutionkernel : :method = evolutionkernel : :N_space; 

will perform the evolution according to eq. @. Other choices are evolutionkernel : :x_space 
and evolutionkernel : : x_space_truncate_in_one_over_L, corresponding to evolution ac- 
cording to eq. ( [TD| ) and eq. (Q), respectively. 

In a similar way, the class alpha_strong has a static data member method, which allows 
one to select the truncation method for a?, e.g. 
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// select truncation method for alpha_strong 

alpha_strong :: method = alpha_strong: :truncate_in_one_over_L; 

will use formula eg .([121), whereas the choice alpha_strong : : exact will use eq. ([TT|). In 
addition the user should set the flavor thresholds and the values of AgcD for the different number 
of flavors, as for example in 

// define values for alpha_strong (in GeV) 
alpha_strong : : charm_threshold = 1.3; 
alpha_strong : : bottom_threshold = 4.7; 

alpha_strong : : lambda_3 = 0.374; 
alpha_strong : : lambda_4 = 0.327; 
alpha_strong : : lambda_5 = 0.226; 

Other classes of general interest provided by the package are the class complex_d correspond- 
ing to complex numbers with double precision as well as the templates vector_template 
<class T> and matrix_template<class T>. The specializations vector_d, vector_c, 
matrix_d and matrix_c define vectors and matrices with real or complex entries in double 
precision. 

4 How to Use the Library 

In this section I give indications how to install and use the program library. Compilation of the 
package will build a (shared) library called libpartonevolution. The user will then link his 
own programs against the library. 

4.1 Installation 

The program library can be obtained from 

http : / / www . f is . unipr . it / ~ stef anw/ part onevolut ion 



After unpacking, the library for the evolution of parton distributions is build by issuing the 
commands 

. / configure 
make 

make install 

There are various options which can be passed to the configure script, an overview can be ob- 
tained with . /configure — help. 

After installation, the shell script partonevolution-conf ig can be used to determine the 
compiler and linker command line options required to compile and link a program with the 
partonevolution library. For example, partonevolution-conf ig — cppf lags will give the 
path to the header files of the library, whereas partonevolution-conf ig — libs prints out 
the flags necesarry to link a program against the library. 
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4.2 Writing programs using the library 



Once the library is build and installed, it is ready to be used. Here is a small example program, 
which defines a parameterization at the scale Qo = 2 GeV and calls the evolution routine for 
(2= 10GeVandx = 0.01 . 

#include <iostream> 

#include "partonevolution/partonevolution . h" 

int mainO 
{ 

using namespace pdf; 

// select truncation methods 

alpha_strong: : method = alpha_strong: :truncate_in_one_over_L; 
evolutionkernel : :method = evolutionkernel : :N_space; 

// define a toy model with 4 flavors 
alpha_strong: : charm_threshold = 0.0; 
alpha_strong: :bottom_threshold = lElO; 
alpha_strong: : lambda_4 = 0.250; 

// define a parton distribution 
double QO = 2.0; 
double A_u[l] = { 35/16.}; 

double alpha_u[l] = { 0.5 }; 
double beta_u[l] = { 3.0 }; 
int eta = -1; 

partondistribution f = 

partondistribution (1, QO, A_u, alpha_u, beta_u, eta) ; 

// evolution 
double Q = 10.0; 
double X = 0.01; 

double res = evolve_nonsinglet (f , Q, x) ; 

// print out result 

cout << "x*f(x,Q^2) = " << x*res << endl; 

} 

After compilation and linking against the partonevolution library, one obtains an executable, 
which will print out 

x*f(x,Q"2) = 0.247242 

This program defines the Wy-distribution 

xu,{x,Ql) = ^/•5(i-^)3-0 (31) 
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of the toy model of ref. [g9|]. The numerical difference of our value with the value 0.24723 
quoted in ref. is well below the estimated accuracy given in ref. [^. 



4.3 Documentation 

The complete documentation of the program is inserted as comment lines in the source code. 
The documentation can be extracted from the sources with the help of the documentation system 
"doxygen" []30[]. The program "doxygen" is freely available. Issuing in the top-level build 
directory for the partonevolution library the commands 

doxygen Doxyfile 

will create a directory "reference" with the documentation in html and latex format. 



5 Checks and performance 

In this section I compare the results of our program with the reference results of Bliimlein et 



al. [^. The authors of ref. [ [290 define a toy model, consisting of four active flavors with 
^QCD ~ 250MeV. At the input scale Qq = 4GeV they take the initial parton distributions as 
follows: 

xuy =AuX^'^{\ -Jc)^, xdv =AdX^'^{\ -x)"^, 
xS = Asx^^-^{l-x)\ xg= Agx^^-^ ( 1 - ^) ^ 

xc = 0, xc = 0. (32) 

The sea S is taken to be SU (3) /^/avor- symmetric and to carry 15% of the nucleon's momentum 
at the input scale. This, together with the usual flavor- and momentum sum rules, fixes the 
constants A^, Aj, As and Ag. The sea is related to the quark singlet distribution 

s= E (q+q) (33) 

q=u,d,... 

by 

S = L-Uy-dy. (34) 

They provide a table with numerical results for the evolved parton distributions at the scale 
— 100 GeV^. They have estimated the numerical accuracy of their results to be 0.02%. I 
have compared the program with the contents of table 1 of [^. Using a 5-point quadrature 
formula the results of the program agree with the numbers presented there within the indicated 
estimated accuracy (0.02%), with some exceptions for the charm distribution and the singlet 
case for small values of x. Using a 10-point quadrature formula for those cases will lead to 
results within the 0.02% uncertainty margin. I have compared both the x-space solution (upper 
part of the table 1 in ref. [P9[[, corresponding to solutions of eq. ( [T^ ) ) as well as the N-space 
solution (lower part of table 1 in ref. [|2^], corresponding to solutions of eq. (|^) in combination 
with eq. ( |T^ for the strong coupling). 

Finally, I would like to give some indication on the required CPU time. The evolution of a 
non-singlet parton distribution with a 5-point Gauss-Laguerre quadrature formula takes about 



12 



1.3- lO^-'* on a standard PC (Athlon, 1.6 GHz). The singlet evolution requires the evalua- 
tion of four anomalous dimensions. In addition there is some overhead for vector and matrix 
arithmetic. The singlet evolution is therefore more expensive in terms of CPU time. The evo- 
lution of the singlet combination takes about 6 ■ 10^5, again with a 5-point Gauss-Laguerre 
quadrature formula on the same PC. These numbers apply to the N-space method. Simulating 
the x-space evolution method, the evolution of a non-singlet parton distribution with a 5-point 
Gauss-Laguerre quadrature formula takes about 1.3 ■ 10^5. In the singlet sector, the evolution 
kernel has to be approximated by evolutions over small intervals. Subdividing the evolution 
interval into 1000 subintervals, the evolution of the singlet combination takes about 270 ■ 10 ■^i', 
again with a 5-point Gauss-Laguerre quadrature formula on the same PC. 



6 Summary 

In this paper I described the program library "partonevolution". This library can be used to 
evolve parton distributions given by parameterizations at a scale Qq to the desired scale Q. The 
program uses the Mellin transform method with an optimized contour. Due to this optimized 
contour the program needs only a few evaluations of the anomalous dimensions and is therefore 
extremely fast. The library is therefore suited to be used in analysis involving parton distribu- 
tions with errors, where a functional integration over a set of pdf 's is needed. Furthermore, once 
the three-loop anomalous dimensions are known, the extension of the algorithms to NNLO is 
straightforward. In addition, the program can also be used to simulate x-space evolution pro- 
grams. 



A Evolution operators 



In this appendix the formulae for the evolution operators are collected. I start with the N-space 
method. The evolution operator, which solves eq. @ is given in the non-singlet case by [Tl|] 



^0 



1 + 



2|3o 



(35) 



Here, T] = — 1 corresponds to the combinations q — q, while T] = +1 should be used for the non- 
singlet combinations like {u + U) — {d + d). I have expressed the evolution operator in terms 
of the anomalous dimensions 7^ and y\ , which are related to the coefficients of the evolution 
kernel by Pf = — ^Y^. The explicit expressions for the anomalous dimensions can be found in 
[p. For the singlet case one finds [Q, [TT] ]: 



/ 

- as{Ql)-as{Q'] 
\ 

+ (+--) 



2po 



pZ pZ 



<Ql) 



pzfpz^ 



(36) 
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with 



0,qg >0,gq 



Pi 



(37) 



The evolution operator solving eq. ( [T0| ) (e.g. "x-space method") reads in the non-singlet sector 

OH 



Po + Pia.(e^) 



(38) 



In the singlet case no closed form is known and the evolution operator has to be obtained as a 
sequence of evolutions over a small interval: 



a 



(n-l) (n-2) 



(39) 



where 

siQo) + {(^siQ^) — (^siQo)) j/n. In a small interval the evolution operator for the 
singlet case can be approximated by [|T^ 

(a.(e'),a,(e^)) 



1 



V 



In 



2pi VPo + Pia.(e^ 

2^1 1,1 + 



2|3o 



,2 + 



2po + ?t^;-?ti 



2|3o ' Po 



-a.(eo)2i^i 1,1 + 



2|3, 



,2 + 







2P0 PO 



+ ( + 



(40) 



Finally, the evolution operator for eq. ( [T^ ) (e.g. "x-space and truncated in 1/L") reads for the 
non-singlet case [ |T2| ] 



-) exp 



pi /l-f21nLo l+21nL\ 



Wo 



L2 
^0 



L2 



/2 + 61nLo + 91n2Lo 2 + 61nL + 91n2L\ , 



L3 
-^0 



L3 



yj(ri) 



(41) 



where L, = In2?/A2. In the singlet case no closed form is known and the evolution operator is 
obtained as a sequence of evolutions 

E'{Q\Ql) = ^^(4o'2(n-i))^'(4-i)'2(.-2))---^'(e(i),e(o)). (42) 
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where = Qo+ {Q^ — Qq) j/n. In a small interval the evolution operator for the singlet 
is then approximated by [jT^ 



case 



pi /l+21nLo l+21nL\ 



L2 



Pf /2 + 61nLo + 91n2Lo 2 + 61nL + 91n2L 



54p6 V 



L3 
-^0 



L3 



'po(2po+^^;-xi) Ilo lU; 



2Pi 



InLo 



InL /Lo 

"p3(4po + ^^+-Xi) ll2--7Tl^7: 



2pi 



1 



1 /^o 



"p2(4po+^^;-xi) Il2 l2Vl 



p5 (6po + X^+-Xi) 
2pf 

p4(6po+^^;-xi) 
if 



8=' 



2P? 



^0 

InLp 

^0 

1 1 



In^L /Lo 
"l3~1 L" 



InL /Lo 
L3" 

-0 ■ 



"p3(6po + ^^+-Xi) \Ll 0\L 



+ (+--) 



(43) 



where 5^ = (X^ — ?ii)/2Po. Note that there is a typo in eq. (8.27) of ref. [|T2[]. The prefactor 

/R,, -X-\~ /R„ 



should read 
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